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IMPLEMENTATION OF LAMINATE THEORY INTO STRAIN RATE 
DEPENDENT MICROMECHANICS ANALYSIS OF 
POLYMER MATRIX COMPOSITES 


Robert K. Goldberg 

National Aeronautics and Space Administration 
Glenn Research Center 
Cleveland, Ohio 44135 


SUMMARY 

A research program is in progress to develop strain rate dependent deformation and failure models for the analy- 
sis of polymer matrix composites subject to impact loads. Previously, strain rate dependent inelastic constitutive 
equations developed to model the polymer matrix were implemented into a mechanics of materials based micro- 
mechanics method. In the current work, the computation of the effective inelastic strain in the micromechanics 
model was modified to fully incorporate the Poisson effect. The micromechanics equations were also combined with 
classical laminate theory to enable the analysis of symmetric multilayered laminates subject to in-plane loading. A 
quasi-incremental trapezoidal integration method was implemented to integrate the constitutive equations within the 
laminate theory. Verification studies were conducted using an AS4/PEEK composite using a variety of laminate 
configurations and strain rates. The predicted results compared well with experimentally obtained values. 


LIST OF SYMBOLS 


A- laminate stiffness matrix components 

b vector of total and inelastic strains in system of equations in micromechanics 
effective stiffness matrix components for composite ply in micromechanics 
D inelastic material constant representing maximum inelastic strain rate 
E elastic modulus of material 

e collected inelastic strain terms in solving for effective inelastic strains 

G shear modulus of material 

h k thickness of ply k in laminate 

J \ second invariant of deviatoric stress tensor 

kj fiber volume ratio of composite 

A J.j total effective force resultant components for laminate 

A^y components of force resultants due to inelastic strains for laminate 

N { number of plies in laminate 

n inelastic material constant representing rate dependence of material 
p vector of unknown subregion stresses in system of equations in micromechanics 
Q.j plane stress stiffness matrix components for composite ply in material axes 

Qij plane stress stiffness matrix components for composite ply in structural axes 
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q inelastic material constant representing hardening rate of material 
R matrix of coefficients of unknown stresses in micromechanics equations 
r inverse of R matrix in micromechanics 

Sy compliance matrix components 
Sjj deviatoric stress components 

t current time 

A t time increment 

T v T 2 coordinate transformation matrices for stress and strain 
Z (} material constant representing initial isotropic hardness of material 
ot scaling factor for shear components of K-, effective stress 

P material constanrused in scaling shear components of K 0 effective stress 

E ij strain tensor components 

£ J jj inelastic strain components 

effective inelastic strain 
£ ( ’j midplane strain components for laminate 

engineering shear strain components 
Y // midplane engineering shear strain components for laminate 
0 fiber orientation angle for each ply in laminate 

v Poisson's ratio of material 

Qy back stress component 

inelastic material constant representing value of back stress at saturation 
Ojj stress tensor components 

C w mean or hydrostatic stress 

• quantities with dots above them represent rates 

Subscripts 

Af bottom left subregion of composite unit cell (fiber material) 

Ml bottom right subregion of composite unit cell (matrix material) 

M2 top left subregion of composite unit cell (matrix material) 

M3 top right subregion of composite unit cell (matrix material) 

R1 region of composite unit cell consisting of subregions Af and Am 

R2 region of composite unit cell consisting of subregions B1 and B2 

C 1 region of composite unit cell consisting of subregions Af and B 1 
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C2 region of composite unit cell consisting of subregions Am and B2 
f fiber related material property 

m matrix related material property 

12 in-plane shear stress or strain components in material axis system 
1 1 ,22 normal in-plane stress or strain components in material axis system 
33 normal out-of-plane stress or strain components in material axis system 
x,y normal in-plane stress or strain components in structural axis system 
xy in-plane shear stress or strain components in structural axis system 


INTRODUCTION 

NASA Glenn Research Center has an ongoing research program to investigate the feasibility of developing jet 
engine fan containment systems composed of polymer matrix composite materials. To design such a system, the 
ability to correctly predict the deformation and failure behavior of the composite under high rate loading conditions 
is required. Specifically, the analysis method must have the ability to account for any strain rate dependence and 
nonlinearities that might be present in the deformation response. 

In previous research (ref. 1 ), an inelastic constitutive model was developed to predict the nonlinear, rate depen- 
dent deformation response of ductile, crystalline polymers. The equations were then implemented into a mechanics 
of materials based micromechanics model to enable the prediction of the nonlinear, rate dependent deformation 
response of carbon fiber reinforced polymer matrix composites. The deformation response of several representative 
uniaxial composites was analyzed for various fiber orientation angles and strain rates. The predicted values com- 
pared well to experimental results, and the rate dependence and nonlinear deformation response of the composites 
were captured. 

In this study, classical lamination theory (refs. 2 to 4) has been implemented into the composite micromechanics 
to allow for the analysis of symmetric thin laminates subject to in-plane loading. As part of the development of the 
laminate theory, the computation of the effective ply level inelastic strains has been modified to fully incorporate the 
Poisson effects into the calculations. Verification studies have been conducted using a representative composite 
composed of AS-4 carbon fibers embedded within a PEEK (polyetheretherketone) thermoplastic matrix. 

In this paper, first a review of the constitutive model used to compute the inelastic deformation response of the 
polymer is presented. Next, the micromechanics model used to compute the effective properties and response of the 
composite plies is discussed, including modifications made in the computation of the effective ply level inelastic 
strains. Finally, the theoretical development and numerical implementation of the incorporation of lamination theory 
into the micromechanics is given, and verification studies are discussed. 


POLYMER CONSTITUTIVE MODEL 
Overview 

The Ramaswamy-Stouffer state variable model (ref. 5), which was originally developed to analyze the visco- 
plastic deformation of metals above one-half of the melting temperature, has been modified to analyze the rate 
dependent, nonlinear deformation response of ductile, crystalline polymers. There is some physical motivation in 
utilizing state variable models that were developed for viscoplastic metals to analyze the nonlinear deformation 
response of ductile polymers. For example. Ward (ref. 6) defined the “yield stress” in polymers identically to how 
researchers (ref. 5) have defined the “saturation stress” in metals. In both cases, the respective terms were used to 
define the stress level where the stress-strain curve in a uniaxial tensile test became flat and the inelastic strain rate 
equaled the total strain rate. Furthermore, Ward (ref. 6) defined an “internal stress’ in polymers to represent the 
orientation dependent resistance to molecular flow during inelastic straining. This concept is similar to the back 
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stress ' in metals, which is defined as the orientation dependent resistance to slip resulting from the piling up of dis- 
locations under a shear stress at a barrier (ref. 5). 

In state variable constitutive equations, a single unified strain variable is defined to represent all inelastic strains 
(ref. 5). Furthermore, in the state variable approach inelastic strains are assumed to be present at all values of stress; 
there is no defined yield stress. State variables, which evolve with stress and inelastic strain, are defined to represent 
the average effects of the deformation mechanisms. 

Several limitations and assumptions have been specified in the development of the constitutive equations. Small 
strain conditions are assumed and temperature effects are neglected. The nonlinear strain recovery observed in 
polymers on unloading is not simulated, and phenomena such as creep, relaxation and high cycle fatigue are not 
accounted for in the equations. The equations are likely to be only valid for ductile, crystalline polymers. 


Flow and Evolution Equations 


In the modified Ramaswamy-Stouffer model, the inelastic strain rate, efj, is defined as a function of the 
deviatoric stress, and tensorial internal stress state variable Q - in the form; 


efj = D ( ,exp 


1 

z 2 1 

n 

S ti 

2 

(3K 2 j 


^ J 


( 1 ) 


where D (y Z () and n are material constants and K 2 is defined as follows: 


2 ( S V Qj ) ( S ij Qj ) 

The elastic components of strain are added to the inelastic strain to obtain the total strain. The following relation 
defines the internal stress variable rate: 


2 ./ , 

— ~ (3) 

w here q is a material constant. £l m is a material constant that represents the maximum value of the internal stress, 
and is the effective inelastic strain. The internal stress is assumed to be equal to zero w'hen the material is in its 
virgin state. The values of the material constants are determined in the manner discussed by Stouffer and Dame 
(ref. 5) and Goldberg (ref. 1 ). 

The hydrostatic stress state has been found to have a significant effect on the yield behavior of a polymer (ref. 7). 
Bordonaro (ref. 8) indicated a possible w ay of accounting for such effects in a state variable constitutive model was 
to modify the effective stress terms. In this work, pressure dependence is included by multiplying the shear terms in 
the K 2 invariant in equation (2) by the following correction factor: 


a = 




(4) 


In this term, c m is the mean stress, J 2 is the second invariant of the deviatoric stress tensor, and (3 is a rate indepen- 
dent material constant. Since only uniaxial data were available for the polymers which are considered in this study, 
the value of the parameter (3 was determined empirically by fitting data from uniaxial composites with shear domi- 
nated fiber orientation angles, such as [15°]. Future efforts will concentrate on developing more systematic methods 
of determining the value of (3. 
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COMPOSITE MICROMECHANICAL MODEL 


Overview 

Micromechanics techniques are used to predict the effective properties and deformation response of the indi- 
vidual plies of the composite laminate. As will be described later, laminate theory is then used to compute the effec- 
tive response of the entire composite laminate. In the micromechanics, the effective properties and response are 
computed based on the properties of the individual constituents. The common procedure is to analyze a unit cell of 
the composite, the smallest material unit for which the response can be considered representative of the response of 
the overall composite. 

The composites are assumed to have a periodic, square fiber packing with a perfect interfacial bond. Small strain 
conditions are assumed, and temperature effects are once again neglected. The fibers are assumed to be transversely 
isotropic and linearly elastic. The matrix is assumed to be isotropic, with a rate dependent, nonlinear deformation 
response computed using the constitutive equations described in the previous section. 

The micromechanics method developed in this study is similar to and based on the approach originally proposed 
by Sun and Chen (ref. 9) and extended by Robertson and Mall (ref. 10). The methodology is also similar to that de- 
rived by Pindera and Bednarcyk in their reformulation of the Generalized Method of Cells (ref. 11). The composite 
unit cell is defined as a single square fiber (with an area equal to the actual circular fiber) and its surrounding matrix. 
Due to symmetry, only one-quarter of the unit cell needs to be analyzed, as discussed in references 9 and 10. Note 
that in reference 1 1 the entire unit cell is used for the analysis. However, for the square fiber arrays assumed in this 
study, the equations resulting from using the entire unit cell would come out identically to those presented here. The 
relationship between the cross-section of the composite, the unit cell, and the section of the unit cell that is analyzed 
is shown in figure 1. The portion of the unit cell used for the analysis is broken up into four rectangular subregions, 
as shown in figure 2. In figure 2, subregion “Af’ represents the fiber and subregions “Ml,” “M2” and “M3” are 
composed of matrix material. The bottom layer of subregions (“Af ’ and “Ml”) is referred to as Row 1 (R1 ), and the 
top layer of subregions (“M2" and “M3”) is referred to as Row 2 (R2). Likewise, Column 1 (Cl ) is defined as con- 
sisting of subregions “Af and “M2”, and Column 2 (C2) is defined as consisting of subregions “Ml” and "M3.” In 
the material axis system shown in the figure, the “1” coordinate direction is along the fiber direction, the “2” coordi- 
nate direction is perpendicular to the fiber in the plane of the composite, and the 3 coordinate direction is perpen- 
dicular to the fiber in the out of plane direction. 


Derivation of Micromechanics Equations 

The transversely isotropic compliance matrix is used to relate the strains to the stresses, using the following rela- 
tions. Note that in these equations represents the components of the compliance matrix, not the components ol the 
deviatoric stress tensor .v, as in the previous section. 
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Y 12 = 5 66°12 +2£l2 (6 

where the symbol “c, ” represents stress tensor components, the symbol “e, ; ” represents strain tensor components, 
and the symbol "y ” represents engineering shear strain components, all assigned in a Cartesian frame of reference. 
A superscript “I” is used to denote inelastic strains. The subscripts “1 1”, “22” and “33" are used to denote stresses, 
and strains along the coordinate axes, while the subscript “ 1 2” is used to denote in-plane shear stresses and strains. 
Out-of-plane shear stresses and strains are neglected in the current analyses. However, these effects could be incor- 
porated into the analyses using similar methods to those presented above. Out-of-plane normal stresses are included 
in the analysis due to the fact that deviatoric stresses are used in the polymer constitutive equations. If these stresses 
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are not included in the analysis, inaccurate results in the inelastic region are obtained. Plane stress can be specified 
on the global level, but the individual subregions must be allowed to have nonzero stresses in all three normal 
directions. 

The addition of the inelastic strain components to the standard transversely isotropic elastic constitutive law 
allows the incorporation of inelasticity into the constitutive relations. For the fiber, which is assumed to be linear 
elastic, these components are neglected. For the matrix material, which is assumed to be isotropic, S, 3 is set equal to 
S 12 , and S 22 is set equal to ( . 

Assumptions of uniform stress and uniform strain are made within the portion of the unit cell used for the analy- 
sis. Along the fiber direction ( 1 direction in fig. 2), strains in each subregion are assumed to be uniform, and the 
subregion stresses are combined using volume averaging. The in-plane transverse normal stresses (2 direction) and 
in-plane shear stresses (1-2 direction) in subregions Af and Ml are equal, while the equivalent strains are combined 
using volume averaging. The same is true for subregions M2 and M3. The total in-plane transverse normal strains 
and in-plane shear strains in Rows 1 and 2 are assumed to be equal, and the equivalent stresses are combined using 
volume averaging. Similar techniques are used for the out-of-plane normal stresses and strains (3 direction). The 
only difference is that Columns 1 and 2 are used instead of Rows 1 and 2. 

By applying the constitutive relations (eq. (5)) for each subregion, and by utilizing the appropriate uniform stress 
and uniform strain assumptions, the following expressions are obtained. In these equations, the total effective strains 
for the unit cell as well as the inelastic strains in each subregion are considered to be known, and the subregion 
stresses are assumed to be unknown. 



e 22 —£ii +y* 
*^I 1m 


“V*7 £ 22M2 + ( l -^f)^ 2BL 4lM3 

^1 Ini ' ' l)i i 


022R2+L/M 
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Equations (11) to (14) form a system of equations that can be solved for the required subregion normal stresses 
a 33C1 and o 33C > These values are substituted into equations (7) to (10), which are then applied to 
solve for the remaining subregion normal stresses. The uniform stress and uniform strain assumptions are used to 
compute the effective total stresses in the unit cell. 

By using equation (6) and the appropriate uniform stress and uniform strain assumptions, the following expres- 
sions are obtained, from which the subregion in-plane shear stresses are computed. Once again the total effective 
strains for the unit cell and the subregion inelastic strains are assumed to be known. 


Y 12 


[>/^7^66f + ( ! ^|^f) S 66m ] CT i2Ri +(2)|l -^7je flMl 


(15) 


Y 12 


= 5 66m CT 12R2 +2 -^7 e 12M2 + _ M 2M3 


(16) 


In the above equations, the subscript “f ' is used to denote fiber related properties, and the subscript "m” is used 
to denote matrix related properties. Subscripts “Af,” “Ml,” “M2," “M3, “Rl, “R2,’“C1 and C2 are used to 
denote stresses and strains in the appropriate region or subregion of the unit cell. Stresses and strains with no region 
identifying subscript are assumed to represent the total effective stresses and strains for the unit cell. The symbol 
“kf ' ’ represents the fiber volume ratio of the composite. 


Effective Inelastic Strain Calculations 


To compute the effective inelastic normal strains for the unit cell, the following procedure is followed. 
Equations ( 1 1 ) to ( 14) are placed into the following format: 


{b}=mp} 


(17) 


where {b} represents the total and inelastic strains in each equation, [R] represents the matrix of the coefficients of 
the unknown stresses, and { p } represents the unknown stresses in the following order: 
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Equation ( 17) is then inverted to solve for the unknown stresses, as follows: 


{/’M'-p} (19) 

w here [r] is the inverse of the [R] matrix in equation (17). The computed stresses in the “22” and “33” directions are 
substituted into equations (7) to (10) to determine the subregion stresses in the “11” direction. The subregion 
stresses are substituted into the equations defining the uniform stress and uniform strain assumptions to determine 
the total stresses for the unit cell, resulting in the following expression: 
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where [C] represents the effective stiffness matrix, and {<?} represents the collected inelastic strain terms resulting 
from the calculations. The vector {*>}is brought to the left hand side of the equation and the resulting expression is 
then inverted. By comparing the result to the transversely isotropic constitutive relation given in equation (5), the 
effective inelastic strains for the unit cell are computed using the following relation: 


£ {\ ] 


( 

e \ 

< ^22 
T 

> = [S]< 

e 2 

e 33 




where {e 7 } represents the effective inelastic strains for the unit cell and [S] represents the effective compliance 
matrix. This procedure differs from what was described in reference 1 in that in the method described here the full 
Poisson effect is accounted for in computing the effective inelastic strains. 


The effective inelastic in-plane shear strain is determined in an identical manner to what was described in refer- 
ence 1. By inverting equations (15) and (16) and using the uniform stress and uniform strain assumptions, the fol- 
lowing expression is obtained: 



where 



(23) 


In this equation, G, 2f represents the shear modulus of the fiber and G m represents the shear modulus of the matrix. 


IMPLEMENTATION OF LAMINATE THEORY 
Overview 7 

In previous research (ref. 1), only uniaxial composites at various fiber orientation angles were analyzed using the 
micromechanics method described above. However, in actual practice, polymer matrix composites are constructed 
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using multilayered laminates, where each ply can have a unique fiber orientation. To allow for the analysis of mate- 
rials of this type, the micromechanics method described above has been combined with classical lamination theory 
(refs. 2 to 4). As will be described below, the force resultants due to inelastic strains are computed in a manner simi- 
lar to that used to determine thermal and moisture resultants in the classical theory. 

Several assumptions have been applied to the analytical method developed here. As in the previous analytical 
development, small strain conditions are assumed, and thermal and moisture effects are neglected. The laminates are 
assumed to be thin enough that the plane stress assumptions normally used in laminate theory are valid. For the cur- 
rent study, out-of-plane stresses are neglected. However, in the analysis of impact problems, transverse shear 
stresses could be significant. Therefore, future efforts will involve adding the ability to incorporate transverse shear 
stresses to the laminate theory. To simplify the theoretical development, symmetric laminates under in-plane loading 
conditions are assumed, and bending effects are neglected. In the future, the ability to analyze unsymmetric lami- 
nates will be added to the theory. 


Basic Equations of Laminate Theory 

In the laminate theory for symmetric laminates, the total force resultants {iV} are related to the midplane strains 
{s 6 *} and the force resultants due to inelastic strains {A^} by the following expression. 


X' 



A] 2 
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r'-i 

e x 
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► = 

A ±2 

A22 

^2^ 

e° y - 



_ A 13 

^23 

^33 _ 

Yxy 


(24) 


where the subscripts “x” and "y” represent normal force resultants, stresses and strains along the “11” and “22” 
directions in the structural (laminate) axis system, and the subscript “xy” represents in-plane shear force resultants, 
stresses and strains in the structural axis system. Note that engineering shear strains are used in the analytical 
development. A schematic demonstrating the difference between the material axis system and the structural axis 
system is shown in figure 3. In the figure, a ply at an arbitrary fiber orientation angle 0 is shown. The material axis 
system “1-2” is displayed, where the “1" direction is along the fibers. The global “X-Y" structural axis system is 
also displayed. All of the plies in the laminate are oriented at various angles related to the structural axis system. 

To compute the laminate stiffness matrix [A], the following procedure is used. First, the effective elastic con- 
stants of the composite in the material axis system are computed using the following equations, which were origi- 
nally developed by Murthy and Chamis (ref. 12). These equations have been found to yield more accurate results 
than those obtained by using the effective stiffness and compliance matrices determined in equations (20) and (21). 
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In these equations, “E" represents the elastic modulus, ‘V’ represents the Poisson's ratio, “G” represents the shear 
modulus, kf represents the fiber volume ratio, the subscript *f represents fiber related quantities and the subscript 
“m" represents matrix related quantities. 

The plane stress stiffness matrix [ Q ] for the composite ply in the material axis system is determined next. The 
matrix has the following format: 


[Q] = 


Q ii 

012 


0 


012 0 
022 0 
0 066 


(30) 


where 


0n = 


1 - Vt^V 


1 022 ~ " 
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12 V 21 




1 — v 12 V 21 1 — V 12 V 2I 

The plane stress stiffness matrix for each ply is transformed into the structural axis system by using the relation: 

[C] = [7',r 1 [e][r 2 ] 


(31) 
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(33) 


and m cos(0) and n=sin(0), where 0 is the fiber orientation angle for each ply. Note that the symbol n has a different 
meaning here than in the polymer constitutive model, where the symbol n represents a material constant. The trans- 
formation matrix f Tj ] is also used to transform the stresses from the structural axis system to the material axis sys- 
tem. The transformation matrix [T,] is also used to transform the strains from the structural axis system to the 
material axis system. 

To calculate the terms in the laminate stiffness matrix [A (> ], the following summation over all of the plies in the 
laminate is carried out. Note that the summation is carried out in the structural axis system. 

N, _ 

['% ] = X [Cf / f hk (34) 

k = 1 


where N t represents the total number of plies in the laminate and h k represents the thickness of each ply. A similar 
summation in the structural axis system is used to compute the terms in the force resultant [A 7 ] : 


m-Spi 

k~\ 


where {e 7 }^ represents the effective inelastic strain vector in the structural axis system for ply k. 


e.n 


hi. 


( 35 ) 
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Incremental Integration of Constitutive Equations 


A quasi-incremental algorithm based on the Heun explicit trapezoidal integration rule (ref. 13) has been devel- 
oped to integrate the polymer constitutive equations. A Runge-Kutta (ref. 13) integration method based on total 
stresses and strains was used in the previous analyses described in reference 1. However, for this work an incremen- 
tal formulation allows for greater flexibility in the application of boundary conditions and simplifies the implemen- 
tation of the laminate theory. 

The following procedure is used within each time step t+At, where t is the time in the previous time step and 
At is the current time increment. First, intermediate values of the inelastic strain rate (denoted by a superscript 1 ) 
in each subregion of the unit cell for each ply of the laminate are computed using equation ( 1 ), where the values of 
the deviatoric stress and internal stress computed for time 1 are used. The first estimate of the inelastic strain incre- 
ments {Ae^-J 1 in the material axis system for the time step are then computed as follows: 

(AeJ-(f-t-Ar)) 1 =(4 <'))‘a/ (36) 

Intermediate values of the internal stress rate (again denoted by a superscript “1") for each subregion of the unit cell 
are then determined using equation (3) (using the value of the internal stress computed for time t). Intermediate val- 
ues of the internal stress at time f+Af are computed as follows: 

Qj,(f + Af) = Q (/ (r) + Qj(OA/ < 37) 


These integrations are similar to a forward Euler type of numerical integration rule. 

Once the inelastic strain increments for each subregion are computed, the effective inelastic strain increments in 
the material axis system for each ply of the laminate are determined using equations (17) to (23 ). The strain^incre- 
ments are transformed into the structural axis system, and the effective inelastic strain force resultants, {AN }, are 
determined using equation (35). 

For strain controlled loading, where one or more terms in the midplane total strain increment vector (Ae^} are 
assumed to be known, a partial inversion of equation (24) is used along with the {AN 1 } to determine any unknown 
midplane strain increments. For stress controlled loading, where increments in the total force vector { AN) are 
assumed to be known, equation (24) is applied directly the compute the midplane strain increments for the laminate. 
The in-plane strain increments in the structural axis system in each ply (x, y, and xy) are assumed to equal the mid- 
plane strain increments for the laminate. After converting the local strain increments to the material axis system for 
each ply, the out-of-plane normal strain increments for each ply are determined. By assuming that each ply is in a 
state of plane stress, the effective elastic constants, the in-plane normal (11 and 22) strain increments and the ply 
level inelastic strain increments are used to compute the out-of-plane total strain increment. 

Given the total and inelastic strain increments for each ply of the laminate, the subregion stress increments in the 
material axis system for each ply are computed using equations (7) to (16). Intermediate values of the total subre- 
gion stresses at time f+A/ are computed as follows: 

c\j(t + At) = aij(t) + Ac)j(t) < 3f 


At this point in the analysis, the “predictor ’ step of the numerical integration is complete. The intermediate val- 
ues of stress and internal stress are then used in equations ( 1 ) and (3) to compute new values of the inelastic strain 
rate and internal stress rate (denoted by superscript “2" in the following equations) for each subregion of the unit 
cell. These approximate values of the inelastic strain increment and total internal stress at time t+At are then deter- 
mined using the following equations. These equations represent a trapezoidal Heun (or second order Runge-Kutta) 
type of integration rule. 


Aef;(f + A/) = 


(4w) 1+ (4( ,+A/ ))' 


'-At 


(39) 
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(40) 


Q,(/ + A /) = %(/) 


ty{t)+a.fj(t + toj\ 


At 


An important point to note is that at this point the values computed using equations (39) and (40) could be com- 
pared to those computed using equations (36) and (37) to see if the difference between the values fall within an error 
tolerance. If the values did not fall within such a tolerance, the size of the time step could be reduced and the values 
could be recomputed. Alternatively, the values computed using equations (39) and (40) could be used as the new 
“predictor” values and the procedure described above could be repeated in an “implicit” procedure. For the current 
study, neither of these methods are employed. Instead, very small time steps are used in order to ensure conver- 
gence. In the future, one or both of the methods described above might implemented into the theory. An improved 
algorithm of this type would improve the flexibility and computational efficiency of the method and allow larger 
time steps to be used. 

Using the revised values for the inelastic strain increments computed using equations (39) for each subregion, the 
effective inelastic strain increments for each ply and inelastic strain force resultants for the laminate are recomputed. 
The midplane strain increments for the laminate are then determined, and the total strains at time t+At are calculated 
by using an expression similar to equation (38). The local ply strains and subregion stresses are recalculated, and the 
revised values of the total subregion stresses are computed using equation (38). The total effective stresses for each 
ply are determined using the uniform stress and uniform strain assumptions for the unit cell, and for strain controlled 
loading these stresses are then combined to obtain the total force resultant at time t+At for the laminate. 


Verification Studies 

To verify the laminate theory implementation, a series of analyses have been conducted using a composite com- 
posed of carbon AS-4 fibers in a PEEK thermoplastic matrix. Tensile stress-strain curves were obtained by Weeks 
and Sun (refs. 14 and 15) for a variety of unidirectional and laminated composites with various fiber orientations and 
laminate configurations at strain rates of 1x10“ and 0.1/sec. Only low strain rate data are examined since the PEEK 
material has only been characterized for relatively low strain rates. However, by conducting verification studies 
using this data, the ability of the laminate theory to capture the rate dependent deformation response of a polymer 
matrix composite has been determined. 

The fiber volume fraction of the AS4/PEEK material is 0.62 (a typical value for this material based on represen- 
tative manufacturer information). The elastic properties of the AS-4 fibers (ref. 12) include a longitudinal elastic 
modulus of 214 GPa, transverse and in-plane shear moduli of 14 GPa, a longitudinal Poisson’s ratio of 0.20 and a 
transverse Poisson’s ratio of 0.25. For the PEEK matrix, the elastic modulus is 4000 MPa and the Poisson's ratio is 
0.40. The inelastic material constants were determined using the procedures outlined in reference 1 and have been 
found to be as follows: D () = lx^/sec, n = 0.70, Z 0 =630 MPa, q = 310, Q m = 52 MPa, |}= 0.00. The value of the 
shear correction factor coefficient P is different from what was used in reference 1. This value is empirically deter- 
mined by fitting data from uniaxial composites with shear dominated fiber orientation angles. Since the method of 
computing the effective elastic constants and effective inelastic strain of the unit cell has been revised since the 
analyses discussed in reference 1 were conducted, it is not surprising that the value of the shear correction coeffi- 
cient also changed. 

Experimental and predicted stress strain curves at both strain rates are shown in figures 4 and 5 for unidirectional 
composites with fiber orientation angles of [30°] and [45°], respectively. These analyses were carried out in order to 
ensure that the laminate theory could still predict the deformation response of unidirectional laminates. As can be 
seen in the figures, the results computed using the nonlinear laminate theory compare well qualitatively to the 
experimentally obtained results for both strain rates. In particular, the nonlinearity and rate dependence observed in 
the experimental curves is captured by the analytical predictions. Quantitatively, the stresses in the elastic range are 
slightly underpredicted, which indicates that the transverse Poisson’s ratio or shear modulus used for the fibers 
might not be correct. The stress level at which the stress-strain curve becomes nonlinear is somewhat overpredicted. 
This could be an artifact of the way the polymer constitutive equations are implemented into the micromechanics or 
the way that the effective inelastic strains are computed. Alternatively, there could be some damage in the experi- 
mental specimens which is not captured by the theory. However, the predicted results are still within approximately 
ten percent of the experimental values at all points in the stress-strain curve. Furthermore, the comparison between 
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the experimental and predicted results improves at higher strain levels in the inelastic range. For failure and life pre- 
diction, correctly predicting the stress levels at the end of the stress-strain curve is most critical. 

To examine the ability of the nonlinear theory to predict the behavior of multilayered symmetric laminates, sev- 
eral laminate configurations were analyzed. Experimental and predicted stress strain curves at both strain rates are 
shown in figures 6 and 7 for composites with layups of [+30°]2 S and [+47°]2 S - Qualitatively, once again the rate 
dependence and nonlinearity observed in the experimental stress-strain curves is captured by the analytical predic 
tions. Quantitatively, for the [±47°] 2s laminates the elastic stresses are slightly underpredicted, and the stress level at 
which the curv es become nonlinear is somewhat overpredicted. The causes of these discrepancies are most likely the 
same as discussed above. However, the stress levels at the end of the stress-strain curve are predicted almost exactly, 
which is most critical for failure and life prediction. For the [±30°] 2s laminates, at a strain rate of 0.1/sec the com- 
parison between the experimental and predicted curves follows the same trends similar as discussed above. 

However, for the lower strain rate, the stress level at which the stress- strain curve becomes nonlinear is more 
significantly overpredicted than is the case for the other analyses. Furthermore, the stresses are overpredicted in the 
entire inelastic range, though there is again some improvement in the predictions at the higher strain levels. In this 
case, there is most likely some significant damage in the experimental specimens which is not captured by the 
theory. 

Experimental and predicted curves at a strain rate of 0. 1/sec are shown in figure 8 for a laminate with a layup of 
[±15°] The predictions shown in the figures compare well qualitatively and quantitatively to the experimental 
results. The important point to note in these results is that the experimental stress-strain curve is highly linear due 
to the fiber dominance of the laminate orientation. This behavior is correctly predicted by the analytical model. 

Even with a nonlinear constitutive model for the polymer matrix constituent, for situations where the deformation 
response is linear the stresses in the matrix are correctly predicted to be low enough for the inelastic strains to be 
negligibly small. 


CONCLUSIONS 

In this study, several modifications to a previously developed methodology designed to predict the nonlinear, rate 
dependent deformation response of uniaxial polymer matrix composites have been carried out. First, the mechanics 
of materials based micromechanics equations used to predict the effective response of a composite ply have been 
modified to improve the calculations of the effective inelastic strains. Furthermore, the nonlinear, rate dependent 
constitutive equations used to model the matrix constituent and the micromechanics model have been combined 
with classical laminate theory to allow' for the analysis of multilayered symmetric laminates. Furthermore, a quasi- 
incremental algorithm based on the Heun trapezoidal integration rule has been developed to analyze the time depen- 
dent deformation of the laminates and integrate the polymer constitutive equations. Predicted stress- strain curves for 
an AS4/PEEK composite compare well to experimentally obtained values for a variety of laminate orientations and 
strain rates. 

Future efforts will include expanding the laminate theory formulations to account for unsymmetric laminates and 
laminates subject to bending loads. Furthermore, the ability to include the effects of transverse shear stresses will be 
added to the laminate theory. The analysis method will also be implemented into a transient dynamic finite element 
code to allow' for the simulation of ballistic impact tests on composite structures. High strain rate experiments will 
be conducted on a representative polymer matrix composite, and the deformation model will be characterized and 
validated for high strain rate conditions. Furthermore, the deformation model will be extended into the large defor- 
mation regime, and the developed techniques will be extended to the analysis of woven composites. 
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Figure 1 . — Schematic showing relationship between composite ply, unit cell, and portion of unit cell which is 
analyzed for micromechanics model. 
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kf = fiber volume fraction 

Af = fiber subregion 

Ml, M2, M3 = matrix subregions 
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Figure 2. — Geometry and layout of portion of 
composite unit cell which is analyzed for 
micromechanics model. 
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Figure 5— Model predictions for AS4/PEEK [45°] laminate at strain rates of 1x10“ 5 /sec 
(IxlO -5 ) and 0.1 /sec. 
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Figure 6.— Model predictions for AS4/PEEK [±30°] 2s laminate at strain rates of 1x10“ 5 /sec 
(1x1 0 -5 ) and 0.1 /sec. 
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Figure 7.— Model predictions for AS4/PEEK [±47°] 2s laminate at strain rates of 
1x10 -5 /sec (1 xIO -5 ) and 0.1/sec. 
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